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A Gaussian operator representation for the many body density matrix of fermionic systems, 
developed by Corney and Drummond [Phys. Rev. Lett, 93, 260401 (2004)], is used to derive 
approximate decoupling schemes for their dynamics. In this approach the reduced single electron 
density matrix elements serve as stochastic variables which satisfy an exact Fokker-Planck equation. 
The number of variables scales as ~ N 2 rather than ~ exp(iV) with the basis set size, and the 
ON ■ time dependent Hartree Fock approximation (TDHF) is recovered in the "classical" limit. An 

approximate closed set of equations of motion for the one and two-particle reduced density matrices, 
provides a direct generalization of the TDHF. 

I. INTRODUCTION 

The study of interacting many-body systems remains one of the most active research fields in physicsi^^. The 
main computational challenge in electronic structure calculations is that a basis of many-body states, spanning the 
, Hilbert space, becomes exponentially large with system size^. In many cases, however, one is not interested in all 
the information contained in the full many-body wave- function (or density matrix). Experimental observations are 
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typically related to expectation values of certain one and two-body operators. This allows for various approximation 
schemes, which focus on correlation functions to retain reduced information on the state of the system. This is the 
basis for the density-functional-theory (DFT) and its time-dependent extension (TDDFT)&£»£»2ii£. 

The many body hierarchy of reduced density matrices can be truncated systematically by making a cumulant 
expansion of the wavefunction. Extensive work had been devoted to reduced description in terms of the coupled one 
, and two-body density matrice o 11 ' 12 . The main difficulty in the direct computation of reduced density matrices has 
been the N-representability problem, namely the lack of exact conditions that guarantee that a given reduced density 
O ! matrix can be obtained from an N electron wavefunction. 

In the TDHF method 1 * ' 4 ' 7 ' 8 ' 13 ' 14 ' 15 ' 16 the state of the system is assumed to be given by a single Slater determinant, 
resulting in a closed system of equations for the occupied single-particle orbitals. This is equivalent to truncation at the 
£N) . level of the reduced single-particle density matrbsi 4 -. TDHF is commonly used as a simple, affordable, approximation 
J> ' for the electronic excitations and the optical responce of a system 1 ^. TDDFT response functions have the same 
formal structure and simplicity as TDHF ones, except that the burden of the many body problem is shifted into the 
' construction of the functiona l 17 ' 18 . 

The TDHF method can be recasted as a system of equations for a set of coordinates, p a p — (c^cg) 13 ' 19 , which 
describe the reduced, single particle density matrix. These can be used to calculate expectation values of any single 
body operator, such as the optical polarization. Note that if the orbital space contain n occupied and m unoccupied 
=~ ■ orbitals, only nm out of the (n + m) 2 elements of this reduced density matrix, namely the electron-hole excitations, 
are needed to represent the full density matrix, and calculate the response^ 4 -. The TDHF equations of motion are 
approximate, and systematic extentions are not obvious due to the absence of a simple small parameter. 

Recently, an exact phase space representation for the many-body fermionic density matrix was developed by Corney 
and Drummon d 20 ' 21 ' 22 . This method is similar in spirit to the coherent states in Hilbert space approach of Cahill 
and Glauber—, since both employ an overcomplete basis set. In Corney and Drummond's approach the many body 
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density matrix is expanded as an ensemble of Gaussian operators A(n) [Eq. |[IU])] in the many-electron phase space, 

p= dnP(n,t)k(n), (1) 



where P(n, t) is a classical probability distribution for the matrix elements n a p. This density matrix is assumed to be 
block diagonal in Fock space, so that coherences are only allowed between states with the same number of particles, 
which is adequate for most applications. Coherences between states differing by an even number of particles, resulting 
in anomalous correlations, can be included as wel L 20 ' 21 ' 22 Grassmann variables, which are essential for the coherent 
state representation in Hilbert space, are avoided, thus providing a more intuitive physical interpretation of various 
quantities. The exact time evolution of the probability distribution P(n, t) for fcrmions with two body interactions 
is described by a Fokker-Planck equation. The dynamics of the many-body system is thus mapped onto an ensemble 
of stochastic trajectories n a p(t) in real (or imaginary) time. Gaussian phase space operator representation have been 
also used to derive the Hartree Fock Bogoliubov equations for Bose Einstein condensates — . 
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So far the Gaussian phase-space representation (GPSR) for fermions has mostly been used to study the ground states 
of Hubbard-like model o 25 i 26 i 27 . Exceptions are the study of a mixed boson-fermion model of molecular dissociation 20 
and a non interacting system^ 2 .. The GPSR can be used to study the time evolution of a system simply by using the 
Fokker-Planck equation to propagate P(n, t) in time. Therefore, this method maps the many body dynamics onto 
an evolution in a classical parameter space. The goal of the current paper is to show that the GPSR can be used to 
derive a new approximation scheme, which naturally connect to, and go beyond the TDHF. The possibility of using 
GSPR to solve the sign problem in imaginary time calculations is under debated However, in the real time dynamics 
considered here complex phases are unavoidable. 

The Gaussian parameters n a p play a role similar to the reduced single electron density matrix elements p a p in 
TDHF. P(n) can be viewed as a distribution of the matrix elements of the reduced single electron density matrix. 
We show that there exist a form of the stochastic equations for n a p(t) whose deterministic part, the non random 
drift terms, coincides with the TDHF equations. It should be noted that the Gaussian operator basis is overcomplete, 
and thus allows for many equivalent forms of the stochastic equations. This freedom may be used to simplify the 
implementation of the method. 

A practical, exact, numerical scheme for computing the stochastic trajectories in real time is yet to be developed. 
We use the GPSR to construct new types of approximations for the dynamics of excitations. By assuming that the 
probability distribution of the parameters n a p is Gaussian, we obtain a closed system of equations for the single-particle 
and two-particle reduced density matrices, which is a direct extension of TDHF. All approaches for performing many 
body computations by using single and two particles density matrices suffer from the N representability problem, i.e., 
it is not guaranteed that the approximate reduced density matrices can be derived from an N electron wavefunctionii. 
The GPSR of the density matrix in Fock space implies that N has a distribution. Developing constraints that restrict 
Eq. |(T]) to a pure state described by a wavefunction with N electrons will be of interest for molecular applications. 
This limitation is less severe for large systems where a distribution of N makes a minimal effect, or for open systems 
like molecules in junctions. 

In Sec. [TT] we present the TDHF equations of motion of the reduced density matrix. This sets the stage for the 
other methods. In Sec. IIHI we present the GPSR of Corney and Drummon d 20 ' 21 ' 22 and derive the Fokker-Planck 
equation. In Sec IIVI we derive an expression for the probability distribution of the number of electrons, and examine 
the conditions which ensure that a representation correspond to a state with a given particle number. In Sec. fVl we 
develop an approximate truncated hierarchy, which goes beyond TDHF, by assuming that the probability distribution 
P(n, t) in Eq. (Q} has a Gaussian form. In Sec. I VII we present stochastic equations of motion, which are equivalent to 
the Fokker-Planck equation, and compare them to the TDHF equation derived in Sec. [ill Our results are summarized 
in Sec. ECO 



II. THE TIME DEPENDENT HARTREE-FOCK APPROXIMATION 



We consider a many-fermion system with two-body interactions, whose Hamiltonian is given by 



a/3 aP'yS 



(2) 



The indices a, (3, 7, 5 denote an orthogonal one particle basis of spin orbitals. The creation and annihilation operators 
satisfy the Fermi anti-commutation rule 



C a C/3 + C/3C Q = S a p. 



(3) 



Without restricting the generality, the two-body interaction V a p^s is taken to be anti-symmetric with respect to 
permutation of the indices a and j3 or 7 and 5. 

Our goal is to derive an equation of motion for a reduced single particle density matrix 



Pa0 = (cic () ) = 

where |$(t)) is the many body wavefunction. We start from the Heisenberg equation 

PeC =imt)\ \H,ctcA !$(*)>, 



(4) 



(5) 



where we work in units such that h = 1. 
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The commutator in Eq. ([5]) is easily evaluated, leading to 

7 a/3jS 

- (*(t)| ct £t 6 7 c c |*(t)) «S 64 + (*(*)| ctcfacg |*(t)) <5< Q - (*(t)| cjct ^gj |*(t)) ^} . (6) 



The TDHF approximation assumes that the state of the system is given by a single slater determinant |$(t)), at 
all times. Using Wick's theorem, we replace the two-body matrix elements with products of single particle matrix 
elements, 

($(t)| c^c^CjCs |$(t)) = p a &P0 1 - p ai P0&- (7) 
Substituting Eq. Q into Eq. © results in the TDHF equations for the reduced density matrix, 

Pe£ = —i ^ {tC-yPe-f — t-feP-ft) ~ * PapP^C, (Vaf/3e ~ V la /3e + V 1<xe p ~ V aJe p) 
7 a/3 7 

— I ^ PafiPeS {VaQSP ~ V a Qf38 + Vta05 ~ V^aSp) ■ (8) 
a/38 

Eq. ijHJ) can be further simplified using the antisymmetry of V a p^s to permutations of a and /?, or 7 and 5. This gives 

Pet = —i X] (HlPei ~~ ^7e/°7C) - 4 * Pa/3P~(t Va~ff3e — P a0 P n&V at& p ■ (9) 

7 a/37 Q /3<5 

Eq. @ implies that can be viewed as classical oscillator coordinates, which follow a deterministic trajector y 13 ' 19 . 
The TDHF equations will be generalized in Sec. [V] using the phase-space representation of a fermionic system, which 
retains the same number of variables, TV 2 , as the TDHF but treats them as stochastic coordinates. We must then 
work with their distribution rather than with deterministic trajectories in that space. 



III. GAUSSIAN PHASE-SPACE REPRESENTATION FOR FERMIONS 



In this section we briefly present the main results of the GPSR, without proofs, following Refs. | 2lll22l We start by 
introducing the Gaussian operators 



A(n) = det n : exp (—6* [21 - n T ] c) 



(10) 



where n = I — n is a square matrix of parameters whose size is given by the spin orbitals basis set, and : • • • : denotes 
normal ordering. A normal-ordered product of creation and annihilation operators is one where the creation operators 
are to the left of the annihilation operators. The sign must be changed when two operators are interchanged during 
reordering. (For instance, : c a c^ := —c^c a . ) The operators A have the following useful properties 



Tr 



TrA 

Tr (kclcp 



1, 

n a p, 

nasnpy — n a7 nps, 
n a pA + ^2 n a7 n S p 

7<5 



Ac! 



n a p 



A + n a7 nsp 



7(5 



dA 
dk 



(11) 
(12) 

(13) 
(14) 

(15) 



where Tr denotes the trace is in the many-body Hilbert space. Generally, the Gaussian operators also obey Wick's 
theorem, 



Tr 



x---*wA:l = £(-l) p Tr 



x • • • x Tr 



^^r-l^^r^" 



(16) 
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where a is a vector composed of all creation and annihilation operators, and i>j = lipfj). The sum runs over all the 
distinct pair permutations. It should be noted that the normal ordering convention used in Eq. (|16p is that only the 
operators which do not belong to A are ordered so that all the creation operators are to its left, and the annihilation 
operators to its right. Eq. (|13[) is a special case of Eq. (Till)) . Proofs, as well as several other similar relations, can be 
found in Ref. [2J Equations (fT4"|) and (fT5)) allow to replace the action of creation and annihilation operators on A by 
derivatives with respect to the parameters n a p. 

With the help of Eqs. ([I]) and (Tl2|) we see that the reduced single particle density matrix is given by the first 
moment of the distribution -P(n), 



Pa/3 



J n al3 P(n,t)dn 



(17) 



We introduce the following notation, for a fermion operator A 



A 



j 'dnP(n,i)Tr Akin) = ^TriA(n) 



(18) 



Note that there is a double averaging here. The Tr takes care of the quantum average for a given set of parameters 
n. We then apply a classical average {■ ■ ■) P over the phase space distribution P(n). 

The reduced p-particle density matrix depends on the lowest p moments of -P(n, t). We will use the expressions for 
expectation values of two and three body operators 



ya/3~f,5e<; 



clc* c\c s c e c ( 



{n/3 7 n a s - n ai nps) p , 



(19) 



(20) 



Before turning to derive the equation of motion for P(n, t), it is important to note that the representation Eq. ([T]) 
is not unique. More than one probability distribution P(n, t) may represent the same many-particle density matrix, 
due to the overcompleteness of the phase space representation. This is easily demonstrated by the fact that the many 
body density matrix of a system with M orbitals only depends on the lowest M'th moments of -P(n). The higher 
moments can thus be arbitrarily chosen. In the following we present the simplest derivation of the Fokker-Planck 
equation, which is naturally related to the TDHF. Alternative Fokker-Planck equations, whose solutions constitute 
different representations of the same many body density matrix, are discussed in App. [XJ 

Substitution of Eq. |T]) into the Liouville equation 



dp 
dt 



(21) 



gives 



dP(n, t) 
dt 



A(n)dn 



P(n,t) H, A(n) 



dn. 



(22) 



The identities (fMj) and (fT5|) can now be used to replace the commutator in Eq. 
that end, we decompose the matrix Va^g in the form 



22|) with a differential operator. To 



Va0~f5 ^ Q a~f ,cQ {35,c- 



(23) 



It should be always possible to find such a decomposition. This decomposition and the number of components c are 
not unique, and the latter can be chosen to be arbitrarily large. This decomposition eventually leads to a Fokker- 
Planck equation whose diffusion matrix is manifestly positive definite. For Hamiltonians with two-body interactions 
Eqs pT]) and (| 1 5[) this leads to a differential operator with second order derivatives 



h,a] =^a 



a/3 



a/3 



OA 

dn a f, 



E 

a/37(5 



Vb (1) b (1) 

I , °af3,c D - 



d 2 A 



EE*, 



3. 



d 2 A 



af 3,c 7S,c dna0dn ^ 



(24) 
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Closed expressions for A a /3 and B^ 2 \ starting with the Hamiltonian Eq. ([2]), are obtained by a straightforward 
but tedious calculation. The decomposition (|23ll allows to write the terms with second order derivatives in the form 
of Eq. (23J, with 

= X! Qa/3,cn a <;n £ /3 (25) 

a/3 
a/3 

where fi was defined after Eq. (JTOJ) . Furthermore, a lengthy calculation results in 

A e ( = —i (t^n^ — t le n 1 c i ) — i «a/3«7C (^a7,8e — ^a^e + V^ Qe/ 3 — V^ 7£/ g) 
7 a/37 

- i ^ n a pn eS {VaCSfi - K*c/3<5 + V^ afi s - V^p) ■ (27) 

a/35 

By substituting (pM|) in Eq. (f2"2")l . and integrating by parts, one obtains a Fokker-Planck-type equation for P(n, t). 

(1 2) 

However, one should note that the coefficients A a /3 and B a p are complex valued rather than real, and it is not clear 



at this point that the term with the second order derivatives is positive definite. A representation with real parameters 
can be obtained by using the analyticity of A(n), and treating the real and the imaginary parts of n a p - 

as independent variables. The identity dA/dn a p = dk/dn^j = -idA/d 

n ^ap can then be used to write 



OA _ {x) dA (!/) dA 
Aa ^- A °^ + A «^- (28) 

Here A^ (A^) denotes the real (imaginary) parts of A respectively. A similar relation holds for the terms involving 
second derivatives and for B^ 12 \ 

After separating the real and imaginary parts of A and B, and performing an integration by parts, one obtains 

a/3 \ an af3 ° n al3 / 



a.fl'yS c 



a/3 7 5 c \ on al3 an 'y5 0n al3 0n 1 8 0n al3 0n 1 S ) 



a/3 7 5 c \ c ' n a/3 ' n 75 On a0 On 1 8 On a0 Orl j8 

Eq. (|29p is a positive definite Fokker-Planck equation, and the probability distribution P(n, t) is guaranteed to remain 
non-negative at all times. In the derivation of Eq. (|29[) boundary terms in the integration by parts were neglected. 
When boundary terms are finite, the methods discussed in App. [A] can be used to derive equivalent Fokker-Planck 
equations with vanishing boundary terms. 



IV. STATISTICAL PROPERTIES OF THE NUMBER OF PARTICLES IN THE GAUSSIAN 

PHASE-SPACE REPRESENTATION 

The Hamiltonian conserves the number of electrons. However the distribution function P(n) may represent 
a statistical mixture of states with different numbers of electrons. It is of interest to find conditions which ensure 
that P(n) represent a fixed number a system with N electrons. So far the GPSR was used to study systems where 
knowledge about the distribution of the number of electrons was not crucial. In this section we derive conditions 
which guarantee that P(n) represents a state with a given number of electrons. This will be needed for application 
to isolated molecules. Furthermore, for general states represented by P(n), we derive expressions for the probability 
to find k electrons in the system, P(fc). 
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We start by calculating the distribution of values of N — ^ Q c^Ca, the number operator, 

It will be convenient to compute the generating function 

G(s) = J dNe- lsN f(N) = (e~ isN ' 
With the help of Eq. (TH)) , it is easy to show that 

d 



NA = ^ n aa + ^ n afi r 

\ a OLfJLU 



dn vp 



A. 



As a result 

G(s) = Tr 

By using the Baker-Campbell-Hausdorff formula we then obtain 

G(s) = /exp 



A. 



1 



isA(n) - ^s 2 B(n) 



(30) 



(31) 



(32) 



(33) 



(34) 



where A(n) — tr(n) and B(n) = tr(n) — tr(n 2 ), where tr denotes the trace in the single particle Hilbert space. (Recall 
that a trace in the many body space is denoted by Tr.) 

Assuming that the integrations over n converge, we obtain the probability distribution (I30[) by the inverse Fourier 
transform of Eq. (f3"Tj) 



1 



exp 



{N-A(n))< 
2B(n) 



(35) 



The distribution P(n) will represent a state with a given number of electrons, No, when f(N) = S(N — No). This 
requires that P(n) > only for points n which satisfy 



tr(n) = tr(n 2 ) = N . 



(36) 



Distributions P(n) which do not satisfy this conditions represent statistical mixtures of different number of electrons. 
To calculate the probability to find k electrons, P(k), we note that the Gaussian operators have a block structure 
in Fock space, with coherences only between states with the same number of particles. For instance, the Gaussian 
operator of a system with two orbitals is 21 

A(n) = detn|00) (00| + (n n h 22 +n 12 n 21 ) |10> (10| + (n n n 22 + n 12 n 21 ) |01> (01| +detn|ll> (11| (37) 
+ n 21 |10)(01|+n 12 |01>(10|. 

The probability to find the system with any number of particles can be found by taking the coefficients of the 
corresponding populations and averaging over P(n). This can be done by defining 



Ci = Trp 2 L •••4 ! c QI •■•c ai , 



(38) 



for I = 1, 2, • • • , M, where M denotes the basis size (i.e. the number of orbitals). Only populations of states with at 
least I electrons contribute to Ci . A direct calculation gives 



(39) 
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This relation can be inverted, leading to 

M 



p M=h^-^ k jrw c - (40) 



which is valid for k = 1, 2, • • ■ , M. 

Wick's theorem (TT51) can be used to write 



Ci=l\ (detn [«!,••• ,ai]) P , (41) 

(oi,- ,ai) 

where sum is over all ordered sets of indices (oil, • • • , a;), while n [ai, • • • , ai] is the minor of n obtained by deleting 
all rows and columns except the ones whose indices are ai, • • • , a/. 

Equations (|40|) and (l4Tj) combine to give the probability to find k electrons in the system 

P(fc) = ( M ]T (detnK---,^]^, (42) 

which is valid for k = 1, 2, • • ■ , M. Equation (|42|) is especially simple for the probability that all the orbitals are filled 

P{M) = (detn) p . (43) 
We obtain the probability that there are no electrons (all orbitals are unoccupied) by using the normalization condition 

AI AI 

P(0) = l-5^P(A) = l + 53(-l)' Yl (detn[o 1> »- l a J ]) J , = <detn) p . (44) 

k—l l—l (ai,--- ,oti) 

In this section we have investigated the correspondence between the phase space distribution P(n) and the number 
of particles in the system. This is one aspect of a broader problem, namely, how to constrain P(n) and its phase space 
evolution, so that it will correspond to a density matrix which exhibits a certain physical property. For instance, what 
are the conditions on P(n) so that it corresponds to a pure state, or to a single slater determinant? The identification 
of such constraints is an open question. 

V. TRUNCATING THE HIERARCHY FOR ONE AND TWO-BODY DENSITY MATRICES 

In this section we use the exact GPSR of the many-body dynamics to develop new approximation schemes which 
provide a natural extension of the TDHF approximation. 

The time evolution of expectation values of any operator follows the Heisenberg equation 

^TrpO = iTrp H,6 . (45) 

For Hamiltonians of the form @ , substitution of normally ordered products of equal numbers of creation and annihi- 
lation operators in Eq. (|45p results in the usual many body hierarchy of equations of motion. The iV'th equation in 
the hierarchy gives the time derivative of an ./V-particle reduced density matrix in terms the reduced density matrices 
of up to N + 1 particles. The first two members in this hierarchy are given by 

~Jj.P< = i l X! taf} (P«C^(8e ~ Pef^(a) + V a pjS (MaPrfSse ~ M a f3,8(S le + M ea ,jsS(p ~ M e p^sS( a ) } (46) 

and 

■^■M^u = [Mg^nySfc - M ae ^ v 8 K + Me^vpS^g - M^pSyg] (47) 

a/3 

aPjS 

— 34ca,M7<5<^/3 + yeCP.^jsS^a ~ D4c/3,i/7<5<^a + -MaP^v [Sse^C ~ ^c^c) + M-tQ^S (SfiaSvp — 5 va 8^)\ , 



8 



where p, M and y were denned in Eqs. (IT). (l"9]) and |20|) . 

The hierarchy for the reduced single-particle and two-particle density matrices will be closed by assuming that the 
probability distribution P(n) is Gaussian. Such an approximation is similar in spirit to the Hartree-Fock approxima- 
tion, since it makes an ansatz on the time dependent state of the system and it includes a single Slater determinant 
as a special case. However, it is more general, since it takes some two body correlations into account. Due to the 
lack of a variational principle, it is not guaranteed that the Gaussian approximation in phase space improves upon 
the Hartree-Fock approximation. However, it is likely to do so. 

When -P(n) has a Gaussian form, averages which are cubic in n a p, such as the ones appearing in Eq. (|20p . can be 
expressed in terms of the first and second moments of the distribution. For instance 

{n a Qnp e n lS ) p = p (np t n lS ) p + (np e ) p (n aC n 75 ) p + (n lS ) P {n a ^np e ) p - 2 (n a f) p (np e ) p (rijs) p . (48) 

Substitution of gS]) in Eq. (20]), with the help of Eqs. (H]) and (JXTJ), leads to 

3-a/37,<K — P0eM a ^,5t + Pa6M/3 y ,eC + PpcMat7,e5 + Py5M a p,eC + PaeMpj£5 + PyiM a p,5 e (49) 
+ PpsM al £ t + PacMp-f^Se + P~feM a pX6 ~ 2 {pa(P0eP~/8 ~ Pa8PPeP~f( + PaSPp<;P~/e 
- PaePPCP-yS + PatP-fCPfiS ~ PaCPpSP-ye} ■ 

Equation (I49|) represent one out of many possible decoupling schemes which can be used to truncate the hierarchy 
of iV-body reduced density matrices. It was derived by making an assumption on the form of -P(n). Since it is based 
on an approximate ansatz on the many body density matrix it has an important advantage over other decoupling 
schemes: it ensures that the approximate density matrix is a physically allowed one, and therefore guarantees that 
the decoupling (|49|l will not result in unphysical expectation values of operators. 

We can thus truncate the hierarchy at the p and M. level. Before doing that, we define the deviation of the two 
body correlation functions from its TDHF values 

AM. a p^& = M a (3,~(6 ~ PaSPPy + PayPpS- (50) 

By substitution of Eqs. (17]), (19]) and (49]) in Eqs. (46]) and (47]) we obtain 

PeC = -i'^ik'yPn - tjePyc) — 4j )] Vg-yPePaPP-yC / ] Vq(5pPaPPed-2i / , VgPesAM a p > SC — 2i }^ V^p^AMep^S, 

7 a/37 a 0S a 0S P^S 

(51) 

and 

AMeC/iu = i ^2 iheAM^u + t j( AM el ^ v - t^AMeC^v - tu^AMeZ,^] ( 52 ) 
7 

- 2i VaPe( [PuvPjin - PafiPpu + AM a 0,fiv] + 2i ^ V^s [PcSPCy ~ Pe-yPCS + AM e { n s] 

ap 7<5 

+ 2i ^ V a p 7 Q \2p ajJi AM e p^u + 2pp 1 AM ea ,^ + 2p av AM t p^ 1 + p tl (AM a p,nv + PuvPp^ - pafiPpu)] 

a{3~y 

+ 2% 2J V a p e S [2ppnAM a C,5v + 2p a& AMpQ^ u + 2p l3u AM a C,ti.8 + PCS (AMap.jj.iy + PauPPfi - PafiPPu)] 
aPS 

- 2i ^ y^PtS [2p ( sAM e p ni/ + 2p eS AMpc, lv + 2pp 1 AM^ v + pp v (AM e c n s + PeSPa ~ P^Pcs)] 

P-yS 

- 2i 2J V a u-ys [2pc-,AM ea ,nS + 2p ei AM a c,ij,s + 2p a sAM^^ 1 + p afi (AM e ( n s + PeSPcj ~ P^yPcs)} ■ 

Equations (j5T]) and ((52]) constitute an approximate closed set of equations for the one particle and two particle reduced 
density matrices which extend TDHF equations to include two-body correlations. The TDHF equation (5]) can be 
recovered from Eq. (51"]) by setting AM. a p a s = 0. 

The exact solution of the Fokker-Planck equation ([29]) does not necessarily satisfy our Gaussian ansatz for -P(n), 
which is why Eqs. (|5ip and (52]) are approximate. One consequence is that they do not conserve the number 
of particles. This can be seen by calculating the time derivatives of the moments of the number operator N, as 
evaluated using Eqs. (51"]) and (52]) . We find that the first two moments are indeed conserved 
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However, this is not the case for the third moment 



(54) 



We conclude that the probability distribution f(N), derived from Eqs. (|5Tj) and (|52| . is not conserved. 

It should be emphasized that the approximations leading to the TDHF equations and to the hierarchy derived in 
this section are of a different nature. The TDHF equations are obtained under the assumption that the state of the 
system is a single Slater determinant at all times. In contrast, here we assumed that the system can be represented by 
a Gaussian operator expansion whose probability distribution is Gaussian at all times. This provides a generalization 
of the TDHF method, suggesting that Eqs. (f5"Tj) and (|5"!?)) should lead to more accurate results than the TDHF 
equations. This approximation scheme should be useful for large systems with many particles, which are likely to be 
less sensitive to changes in the distribution of number of particles. 

VI. UNRAVELLING THE FOKKER-PLANCK EQUATION IN TERMS OF ITO STOCHASTIC 

TRAJECTORIES 

A probability density evolving according to a Fokker- Planck equation can be computed numerically by simulating 
an ensemble of trajectories following a stochastic equation of motion 2 ^. This is commonly used to trade off computer 
memory by time. This approach was used to study the ground states of Hubbard-type model o 20 ' 22 i 25 ' 26 i 27 . Here we 
present the form of the stochastic equations of motion that naturally connect with the TDHF equations. 

Stochastic equations of motion which are equivalent to the Fokker-Planck equation ([29"]) are given by 



where we have combined the equations for the real and imaginary parts of n a/ 3 into a single complex equation. 
The noise Wiener increments are real, with Gaussian probability distribution, satisfying (dWc (t)dwj:P (f) \ = 



Scc'SijS (t — t') dt. Equation ((55)) should be integrated using Ito stochastic calculus 3 ^. The many-body dynamics can 
be simulated by repeatedly integrating Eq. ([55]) to create an ensemble of stochastic trajectories. This ensemble, 
together with Eq. (JT7J) , or similar identities, can be used to calculate expectation values. 

The non-random, or drift, component of the stochastic equation of motion (|55[) . given by (|27[) . coincides with the 
right hand side of Eq. ([8]) , as long as we identify n a p with p a p . Such identification can naturally be made when the 
distribution P(n) is narrow, or when the noise part is absent (for instance when there are no two-body interactions). 

It is possible to derive an equivalent form of the stochastic differential equations by using a different stochastic 
calculus^ .. For instance, the Stratonovich calculus results in a similar stochastic equation with A^t™^ ^ A a p- This 
breaks the simple correspondence with TDHF. For completeness, we give an explicit form of the drift terms in the 
Stratonovich scheme in App. [Bj 

The coefficients A a p, and B a g' c , appearing in Eq. ([55|) . depend quadratically on n. As a result, the nonlinear 
stochastic trajectories generated from Eq. (|55p may escape to infinity. This is expected to happen in the vicinity 
of certain phase space regions. This problem may be solved in various ways. When calculating the ground state 
using an imaginary time equation, one can force the stochastic equations of motion to be rea l 20 i 22 ' 25 . This limits the 
dynamics to a small subset of phase space, which is not likely include the most unstable regions. Alternatively some 
combination of the gauge freedoms, described in App. \K\ can be used to obtain an equivalent, more stable, form of 
the stochastic equations. A simple systematic way of obtaining stable stochastic equations for the time evolution of 
general fermionic models is yet to be developed. 



In this paper we have used a Gaussian phase-space representation to generalize a common Hilbert space approxi- 
mation for many-body dynamics, the TDHF, which assumes that the wavefunction of the system is given by a single 
Slater determinant at all times. Here the many body density matrix is represented by an ensemble of Gaussian 
operators whose distribution of parameters satisfies a Fokker-Planck equation. We have shown that the drift terms 
of this Fokker-Planck equation have similar form to the TDHF equations [compare Eqs. (f!?T|) and ©]. This suggests 
that the phase-space representation can serve as a good basis for new approximation schemes, which systematically 
go beyond TDHF. 




(55) 



C 




VII. DISCUSSION 
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The Fokker Planck equation (|29|) maps the original fermion system rigorously into a equivalent set of classical 
oscillators n a p with stochastic dynamics. The number of oscillators is the same as the number of elements of the 
reduced single electron density matrix and it thus scales as ~ N 2 with the basis set size. The TDHF is a mean field 
approximation to this stochastic dynamics whereby the same classical coordinates satisfy a deterministic equation. 
Generalizations of the type given here were conjectured in ref. [19| but the GPSR puts it on a firm theoretical basis. 

Linear and nonlinear response functions e.g. optical spectra have a very different formal structure for quantum 
and classical systems. To calculate the response function to a classical field E(t) one needs to add a coupling term 
E(t)A to the Hamiltonian where A is a dynamical variable, nth order quantum response functions are then given by 
n nested commutations such as ([[A(n), A(t2)] , A(t3)]) which gives 2™ terms (Liouville space pathways). Classical 
response functions have a very different form and require to run groups of several nearby trajectories and carefully 
monitor how they diverge. The responce functions can be expressed in terms of the stability matrices of the stochastic 
trajectories of the free system (without external driving) 31 ' 32 , or by solving the distribution of the driven system and 
combining terms with different orders of interactions (finite field techniques)^ 3 .. An intriguing aspect of the present 
approach is that the quantum response functions are recast rigorously in a classical form. This could result in new 
insights and could suggest new numerical simulation techniques for fcrmions. 

Studying the merits of this new representation will be an interesting future direction. Moreover, the quantum 
response function is one specific combination of the Liouville space pathways. Other combinations represent spon- 
taneous fluctuations and how they are affected by external drivin g 34 ' 35 . Many attempts have been made to connect 
nonlinear response and fluctuations through generalized fluctuation-dissipation relations 36 . The GPSR may shed a 
new light into this issue and provide a classical picture for the connection between quantum response and fluctuations. 

We have constructed an approximate, simple, numerical scheme by further assuming that the probability distribution 
P(n) has a Gaussian form at all times. The hierarchy of equations for reduced many-body density matrices can be 
truncated at the single and two particle level. This yields a coupled set of equations, (|5ip and ([52"]) . which can be 
solved to calculate the time dependence of expectation values of all the one and two-body operators. 

The exact time evolution of the phase-space distribution function, which represents the many body density matrix, 
conserves the distribution of the number of electrons. This is not true for the approximate dynamics of Eqs. (|51[) and 
(|52[) . While the mean number of electrons, and its variance, do not vary with time, higher moments are not conserved. 
This approximation scheme thus suffers from the old N representability problem 12 common to other reduced density 
matrix approaches. Nevertheless, this scheme may be a useful description, most likely for systems with many electrons, 
where expectation values are less sensitive to small variations in the distribution of the number of particles. 

The phase-space representation allows the development of new approximation schemes, by making other assumptions 
on the form of P(n). It may be possible to include constraints that ensure that approximations on -P(n) conserve the 
distribution of number of particles, while still allowing to truncate the hierarchy. More work is needed in order to gain 
understanding on the way in which various assumptions on P(n) manifest thcmeselves on the many body dynamics. 
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APPENDIX A: GAUGE FREEDOM IN THE FOKKER-PLANCK EQUATION 

As is the case for Bosonic coherent states^, the Gaussian operator basis for fermions is overcomplete. This means 
that there exist an infinite number of equivalent forms of the time evolution equations for P(n) which give the same 
expectation values of all physical observables. Below we give a brief description of the various ways of obtaining such 
equivalent representations, following Ref. 22, and examine whether they retain the correspondence between the drift 
terms and the TDHF equations. 

A way of changing the Fokker-Planck equation, without affecting the physical observables, is to add formally 
vanishing operators, such as c^c^c^c-y, to the Hamiltonian, and then use Eqs (TH| and (fT5|) to replace these terms by 
a differential operator. This method, which is only applicable for systems of identical fermions, was termed Fermi 
gauge in Refs. [20II22I A direct calculation shows that such Fermi gauges do not change the form of the drift terms 
Aap in Eq. (|2T))) . but do affect the noise terms. Fermi gauges were used to obtain imaginary-time equations for P 
where all the parameters are real rather than complex valued. 

A different method of obtaining equivalent forms of the equations of motion is termed drift gaug o 22 ' 38 . Here 
one adds a new parameter, which corresponds to a weight given to each stochastic trajectory, thus replacing the 
(uniform) average over trajectories by a weighted average. As a result, a new distribution function is defined in a 
larger parameter space, which includes the weight as a variable, and a Fokker-Planck equation is derived for this new 
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distribution functions. At the same time, the evaluation of averages using the distribution P is replaced by weighted 
averages involving the new weight parameter. It is then possible to modify the drift terms A a j3, while at the same 
time adding suitable noise terms to part of the Fokker-Planck equation related to the new weight parameter, in such 
a way that will not change the values of weighted average s 22 ' 25 . Drift gauges have been used to change the drift terms 
in order to avoid trajectories which run to infinity. 

An equivalent, but different form of the stochastic equations of motion can be obtained by using different forms for 
the decomposition of Q in Eq. (j2"3")l . This will not affect Eq. (f!?9")) . but will change its stochastic trajectory simulation, 
for instance using Eq. (|55[) . However, when combined with drift gauges the choice of solution of Eq. (|23[) may result 
in a different Fokker-Planck equation in the generalized phase space. The freedom to choose different solutions of Eq. 
was termed diffusion gauge in Refs. [2M38I . 

The gauge freedoms may be useful for practical implementation of the GSPR. If a solution of Eq. (j2"9")l leads to a 
probability distribution with long tails, and therefore for the appearance of boundary terms, due to the integration 
by parts, the gauge freedom can be used to generate an equivalent representation whose probability distribution is 
more localized. One can hope to use the many ways to generate equivalent representations, by combining the various 
gauges, to find a representation which does not suffer from boundary terms. 



APPENDIX B: THE STRATONOVICH STOCHASTIC EQUATIONS 

The Stratonovich calculus is an alternative to the Ito calculus used in eq. (|55|) . In this scheme the terms in the 
stochastic differential equation are evaluated at the midpoint of the interval, as opposed to the initial point used in 
the Ito scheme. 

In practice, one solves the finite differences equation 

ntT ] - n % = ^ ^ mld ) dt + E (C> m "V^ (1) + B^ mid )dWP) , (Bl) 

c 

with xi mld = (n( I+1 ) +nW)/2, in order to get n^ +1 ^. This implicit equation can be solved by iteration. It is of interest 
to note that implicit methods tend to show better numerical stability then explicit ones, see Ref. [39 for a review. 

This stochastic integration should lead to the Fokker-Planck equation (|2"9"|) . However, this means that the coefficients 
j^(strat) can na |. k e e q ua j ^ ^he ones appearing in Eq. (I2T[) . This is due to correlations in the noise terms, which 
lead to the second order derivatives with the form ^Jr-B-S-BP, rather than \ -S--3-BBP which would appeared when 

2 on On ' 2 On On 

using the Ito scheme. (We suppressed subscripts for brevity.) 

To compensate for these correlations, the drift terms of the Ito and Stratonovich schemes are related by 

7(5 c ' 75 c ' 

For the Hamiltonian © we find 

j^(strat) _ _■ ^^ n ^ _ t ya n 7f3 ) — i (Vf3SS-/n a7 — VySSgnyp) 
7 7(5 

+ * X/ n 7< 5nQe (V/3je6 + VjffSe) - i X n lS n e p (V e7a s + Vy e Sa) ■ (B3) 
7(5e 7(5e 

We have just shown that the Stratonovich form of the drift term differ from the one obtained using the Ito scheme, 
and therefore, differs from the terms appearing in the TDHF equations of motion. Nevertheless, by construction, 
both stochastic schemes give the same ensemble of trajectories. 
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